The self-assembly and evolution of homomeric protein complexes 
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We introduce a simple "patchy particle" model to study the thermodynamics and dynamics of 
self-assembly of homomeric protein complexes. Our calculations allow us to rationalize recent re- 
sults for dihedral complexes. Namely, why evolution of such complexes naturally takes the system 
into a region of interaction space where (i) the evolutionarily newer interactions are weaker, (ii) 
subcomplexes involving the stronger interactions are observed to be thermodynamically stable on 
destabilization of the protein-protein interactions and (iii) the self-assembly dynamics are hierarchi- 
cal with these same subcomplexes acting as kinetic intermediates. 

PACS numbers: 87.15.km,87.14.ak,81.16.Dn,87.23.Kg 



INTRODUCTION 

A large proportion of proteins are not monomeric in 
vivo, but instead 50-70% of those of known structure ex- 
ist as homomeric protein complexes Q . These com- 
plexes are usually symmetrical with each protein in an 
identical environment This latter constraint limits 
the number of possible symmetries for such complexes. 
Thus, for homomeric tetramers, the two possible sym- 
metries that obey this rule are cyclic (C4) or dihedral 
(1^2) (Fig. [1]). The C4 geometry involves only one type 
of interaction, whereas the D2 complex involves at least 
two self-complementary interactions. Interestingly, dihe- 
dral complexes are over ten times more abundant than 
cyclic complexes with the same number of subunits [3]. 
The origin of this preference seems to be evolutionary, 
namely because self-complementary interactions are eas- 
ier to generate and because the evolution of dihedral 
complexes from a monomer does not have to proceed in 
a single step, e.g. a C2 dimer can be an intermediate on 
the evolutionary pathway to a D2 tetramer. 

Insights into the evolution of homomeric protein com- 
plexes have recently come from a study that compared 
the complexes adopted by homologous proteins [3]. Al- 
though, for most homologues the quaternary structure 
is conserved, those cases where it is not conserved can 
tell us something about the evolutionary relationships 
between complexes of different symmetry. So, for exam- 
ple, it was found that where a tetramer shared an evo- 
lutionary relationship with a dimer, the tetramer always 
had D2 symmetry and in the majority of cases the dimer 
interface was conserved in the tetramer, supporting the 
postulated role of the dimer as an evolutionary interme- 
diate. The study went on to show that the evolutionarily 
older interface was usually larger [3], and so presumably 
had a stronger interaction strength. In addition, mass 
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FIG. 1: (Colour Online) Schmatic depiction of (a) a D2 
tetramer and the possible equilibria involved in its formation, 
and (b) a C4 tetramer. 

spectrometry revealed that when the disassembly of di- 
hedral complexes was induced by changing the solution 
conditions (e.g. through the addition of denaturant) sta- 
ble subcomplexes involving the larger interfaces were de- 
tected in the majority of cases [3]. Thus, the static struc- 
ture of a dihedral complex, in particular the ratios of the 
interface areas, can provide insight into the evolution and 
assembly of the complex. 

Here, we provide a framework for understanding these 
results for dihedral complexes by studying the thermody- 
namics and dynamics of assembly for a simple model of 
the complexes, particularly focussing on the role played 
by the relative strengths of the different interactions. 

As the proteins in oligomeric complexes interact 
through highly specific contact surfaces or "patches" , we 
model the proteins using a patchy particle model used 
previously [G*, 't^ , but with the introduction of an addi- 
tional torsional component to the potential that ensures 
that the patches must not only point at each other to 
interact strongly, but also have the correct relative orien- 
tation. In this model, the repulsion between the particles 
is based upon an isotropic Lennard- Jones potential: 



VLj{r) = 4eref 
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but where the attraction is modulated by an orientational 
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term, V^ng- Thus, the complete potential is 
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FIG. 2: (Colour Online) (a) Free-energy diagram showing the 
dependence of the most stable state of the D2 tetramer system 
on temperature and relative interaction strengths. The lines 
indicate where the equilibrium constants for the association 
reactions in Fig. [TJa) are one, and are solid when the equi- 
librium is between the two most stable forms, and dashed 
otherwise. The diagram is shaded according to whether 
monomers, dimers or tetramers are most stable. The arrow 
indicates a possible evolutionary path from a dimeric to a 
tetrameric complex, (b) Dependence of the final yields of 
monomers, dimers and tetramers from our dynamics simu- 
lations on kT/cAA and esB/^AA- Each pixel represents the 
result of a separate simulation, where each simulation started 
from a random configuration of the 400 particles and was 10^ 
steps long, (c) Equilibrium probability of a particle being 
monomer ic, or in a dimer or tetramer at cbb/^aa =0.2. 
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where Q,i is the orientation of particle i, and 
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where Oaij is the angle between the normal to patch a 
on particle i and the interparticle vector Vij , ^ is a tor- 
sional angle, and the 'max' selects the pair of patches 
that have the strongest interaction for the current geom- 
etry. cTpJ and a^^l are measures of the specificity of the 
patch-patch interactions, and which we here choose to be 
the same for all patches (l5|. By contrast, we allow the 
well-depth of the patch-patch interactions, e^fd to vary 
(eref = max[eQ,/3]). We also assume that interactions be- 
tween pairs of patches not in contact in a complex to be 
zero (e.g. cab = for D2 tetramers). 

To simulate the dynamics of our systems we use the 
'virtual move' Monte Carlo algorithm of Whitelam and 
Geissler [s^ as this generates the diffusional behaviour ex- 
pected of particles and clusters in solution. To determine 
the thermodynamic properties of the system, we analyt- 
ically calculate the partition functions for each state, as 
an ideal gas of clusters with rotational and vibrational 
degrees of freedom [16]. The vibrations are assumed to 
be harmonic and their frequencies are calculated by di- 
agonalization of the Hessian. 

We first consider the case of D2 tetramers, which for 
simplicity we choose to be planar. Fig. [2](a) shows the 
thermodynamics of the system as a function of the ra- 
tios of the interaction strengths cbb/^aa- Note, due 
to the symmetry we need only consider the case where 
^AA > ^BB- At high temperature T (or equivalently low 
^aa) the system is monomeric. At cbb/^aa = 0, as 
the system is cooled, it passes from a 'gas' of monomers 
to a 'gas' of dimers. At non-zero cbb/^aa a transition 
from dimers to tetramers appears, and its transition tem- 
perature increases with 63 b/ ^aa until a critical value of 
^bb/^aa ^ 0.5 is reached beyond which dimers are no 
longer most stable for any value of kT/cAA This 
critical value corresponds to the value of cbb/ ^aa at the 
'triple point' where the three equilibrium lines in Fig. 
[2](a) meet. 

This disappearance of the dimeric state is easy to un- 
derstand. In the monomer to dimer transition the effec- 
tive number of particles decreases by half, and the energy 
decreases by e^A, whereas in the dimer to tetramer tran- 
sition, although the number of clusters again decreases 
by half, the energy decreases by 2eBB- Therefore, the 
value of €bb/^aa for which the monomer to AA dimer 
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FIG. 3: (Colour Online) (a) and (b) Time dependence of the probability of a particle being in a particular state. For (a) 
esB/eAA = 0.6 and kT /caa = 0.05 and (b) CBB/eAA = 1 and kT /caa = 0.04. (c) Dependence of /i, the hierarchicity index, on 
temperature and relative interaction strengths. The state points corresponding to (a) and (b) are marked by a white dot. 



and the AA dimer to tetramer transitions occur at the 
same temperature is expected to be less than one, and is 
predicted to be 0.5, in agreement with Fig. [2ja), if the 
entropy change for both of these 'dimerizations' is simply 
assumed to be equal. 

The reliability of our analytical thermodynamic model 
is confirmed by comparisons to results obtained using the 
Wang-Landau algorithm Q and to our dynamics simula- 
tions. The kinetic yields of monomer, dimer and tetramer 
presented in Fig. EJ^b) mirror the form of the free en- 
ergy diagram (Fig. [2]^a)), showing that near equilibrium 
is reached in the simulations except at very low tempera- 
tures where tetramers fail to form properly. In this latter 
case, the rate of cluster formation is so fast and the rate 
of cluster breakup so slow that all the monomers are used 
up before complete tetramers can form. Instead, many 
particles get trapped in trimers (e.g. Fig. [31(b)). 

We can use Fig. [Sfa) to help us understand the evo- 
lution of a D2 tetramer. The arrow in this figure shows 
a possible evolutionary pathways that takes the system 
from a dimeric state (i.e. where cbb/^aa = 0) to a region 
of interaction space where tetramers are stable. The im- 
portant thing to note is that the result of this evolution 
is a tetramer where the evolutionary newer BB interac- 
tions are weaker. This conclusion is in agreement with 
Levy et al.^s observation that newer patches have smaller 
areas [3] , assuming that the interface area is a reasonable 
proxy for estimating the interaction strength. 

We can also use Fig.[2](a) to consider the effect of addi- 
tives that cause disassembly of protein complexes. If such 
additives destabilize all interactions equally, this process 
would correspond to a vertical pathway in Fig. EJa). 
Therefore, if a complex has a value of cbb/^aa ^0-5, 
a thermodynamically stable dimeric phase will be seen 
for some concentration of the denaturant. The equilib- 
rium probabilities of being in the different states along 
one such pathway is illustrated in Fig. EJc). Again, as 
evolution will naturally take the system into a region to 
the left of the 'triple point', our results are in agreement 
with the mass spectroscopic results of Ref. i3i that de- 



tected subcomplexes involving the stronger interactions, 
and with more traditional studies of tetramer formation 
For example. Fig. EJc) is very similar to results 
presented for phosphofructokinase [l2[. 

Fig.[2](a) also allows us to think about the kinetic mech- 
anisms of tetramer formation. In particular, the region in 
which tetramers are stable can be divided up into three 
subregions, depending on the stability of A A and BB 
dimers. In region N neither dimers are stable with re- 
spect to monomers, and so all intermediates are unsta- 
ble and there will be a nucleation free energy barrier to 
tetramer formation. In region AA dimers, but not 
BB dimers, are stable with respect to monomers, and 
so AA dimers will act as kinetic intermediates. In this 
region, we therefore expect a hierarchical self-assembly 
mechanism to dominate, in which AA dimers first form, 
and which in turn dimerize to form tetramers, rather 
than a mechanism which proceeds by sequential addi- 
tion of monomers. Finally, in region D both A A and 
BB dimers are stable with respect to monomers and so 
all pathways for tetramer formation are downhill in free 
energy. 

To test these predictions, we first look at the time 
dependence of the probabilities of being in the differ- 
ent states. In region the sequential passage from 
monomers to dimers to trimers to tetramers is evident 
(Fig. [31(a)). However, the initial rapid rise in the number 
of tetramers slows down as monomers become depleted, 
and the further formation of tetramers from trimers is 
dependent on cluster breakup (a relatively slow process) 
releasing additional monomers. By contrast, in region 
AA dimers are clear intermediates and there is a steady 
growth of tetramers indicative of formation by dimer- 
dimer addition (Fig. [31(b)). Indeed, this plot has a very 
similar form to results for experimental studies on the 
rate of tetramer formation, such as for phosphoglycerate 
mutase [13] and lactate dehydrogenase [14]. 

We have further analysed how hierarchical the dynam- 
ics are by introducing a 'hierarchicity' index h that we 
define as the fraction of tetramer forming events that 



4 




^ ' ~ ''"i 

0.2 0.4 0.6 0.8 1.0 0.8 0.6 0.4 0.2 



plexes where the newer patches are weaker and where the 
assembly is hierarchical. 

In this paper we have shown how an analysis of the 
dependence of the thermodynamics and kinetics of ho- 
momeric protein complexes on the relative interaction 
strengths can provide a framework for understanding how 
their properties are constrained by their evolution, in 
particular their asymmetry in their interface areas and 
their hierarchical assembly. Although we have focussed 
on tetrameric complexes, the approach is easily generaliz- 
able to larger complexes and leads to similar conclusions. 
Thus, it would be very interesting if diagrams similar to 
the free energy diagrams of Figs.[2fa) andjU^b) could be 
mapped out experimentally by studying in detail how the 
thermodynamics of disassembly depends on the relative 
interface areas in a variety of protein complexes. 

The authors are grateful for financial support from the 
EPSRC and the Royal Society. 



FIG. 4: (Colour Online) (a) Possible equilibria involved in the 
formation of a D4 octamer. (b) Free energy diagram show- 
ing the dependence of the most stable state on temperature 
and interaction strengths, superimposed on the final yields of 
monomers, dimers, tetramers and octamers in our dynamics 
simulations (800 particles, 10^ steps). Note the change in the 
ordinate and abscissa in each half of the diagram. 

occurred by dimer-dimer addition weighted by the frac- 
tional yield of tetramers. It can be clearly seen that the 
region in Fig. [31(c) where h is high corresponds well to 
region H in Fig. [2](a). 

The approach we have put forward in this paper for un- 
derstanding the formation of tetrameric complexes can be 
equally applied to complexes with other numbers of sub- 
units. To illustrate this, we show the free energy diagram 
for the formation of a D4 octamer in Fig. [Ub) calculated 
by a simpler version of the theory used for Fig. [2](a), su- 
perimposed on the yields of different clusters obtained 
from our dynamics simulations. In our model for this 
system, the particles have three patches, two types of in- 
teraction (AB and CC) and for simplicity the octamer 
has a cubic shape (Fig. [H^a)). There are two possible 
hierarchical pathways for octamer assembly either via a 
C4 tetrameric intermediate that is stabilized by the AB 
interactions (region Ht) or via dimers stabilized by the 
CC interactions (region Hjj)^ and this expectation is con- 
firmed by analyses of the dynamics in these regions. In 
comparison to the results for the tetramer, it is notice- 
able that the low temperature region in which incomplete 
assembly leads to poor yields extends to higher tempera- 
ture; this is a general trend for the self-assembly of more 
complex targets. The arrows in Fig. [IJb) illustrate po- 
tential evolutionary pathways of an octamer from a cyclic 
tetramer or a dimer, and which will again lead to com- 
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